# Replication Code For Figure 4 and Results

min <- read.csv("./min.csv")

library(dplyr)
min <- min %>% filter(treat >2) %>% mutate(treat = treat - 3)

#find all of the effect sizes, pvalues, and se
mean4 <- summary(lm(data = min, support ~ treat))$coefficients[2,1]
se4 <- summary(lm(data = min, support ~ treat))$coefficients[2,2]
upper95_4 <- mean4 + se4*1.96
lower95_4 <- mean4 - se4*1.96
upper90_4 <- mean4 + se4*1.645
lower90_4 <- mean4 - se4*1.645
p4 <- summary(lm(data = min, support ~ treat))$coefficients[2,4]

min <- read.csv("./min.csv")

min <- min %>% filter(treat < 3) %>% mutate(treat = treat - 1)
mean3 <- summary(lm(data = min, support ~ treat))$coefficients[2,1]
se3 <- summary(lm(data = min, support ~ treat))$coefficients[2,2]
upper95_3 <- mean3 + se3*1.96
lower95_3 <- mean3 - se3*1.96
upper90_3 <- mean3 + se3*1.645
lower90_3 <- mean3 - se3*1.645
p3 <- summary(lm(data = min, support ~ treat))$coefficients[2,4]

txdata <- read.csv("./tex.csv")
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

mean1 <- summary(lm(data = tx_pos, support ~ treatment))$coefficients[2,1]
se1 <- summary(lm(data = tx_pos, support ~ treatment))$coefficients[2,2]
upper95_1 <- mean1 + se1*1.96
lower95_1 <- mean1 - se1*1.96
upper90_1 <- mean1 + se1*1.645
lower90_1 <- mean1 - se1*1.645
p1 <- summary(lm(data = tx_pos, support ~ treatment))$coefficients[2,4]

mean2 <- summary(lm(data = tx_neg, support ~ treatment))$coefficients[2,1]
se2 <- summary(lm(data = tx_neg, support ~ treatment))$coefficients[2,2]
upper95_2 <- mean2 + se2*1.96
lower95_2 <- mean2 - se2*1.96
upper90_2 <- mean2 + se2*1.645
lower90_2 <- mean2 - se2*1.645
p2 <- summary(lm(data = tx_neg, support ~ treatment))$coefficients[2,4]

min <- read.csv("./min.csv")

min <- min %>% filter(treat < 3) %>% mutate(treat = treat - 1)

mean7 <- summary(lm(data = min, min_ben ~ treat))$coefficients[2,1]
se7 <- summary(lm(data = min, min_ben ~ treat))$coefficients[2,2]
upper95_7 <- mean7 + se7*1.96
lower95_7 <- mean7 - se7*1.96
upper90_7 <- mean7 + se7*1.645
lower90_7 <- mean7 - se7*1.645
p7 <- summary(lm(data = min, min_ben ~ treat))$coefficients[2,4]

min <- read.csv("./min.csv")

min <- min %>% filter(treat > 2) %>% mutate(treat = treat - 3)

mean8 <- summary(lm(data = min, min_ben ~ treat))$coefficients[2,1]
mean8
se8 <- summary(lm(data = min, min_ben ~ treat))$coefficients[2,2]
upper95_8 <- mean8 + se8*1.96
lower95_8 <- mean8 - se8*1.96
upper90_8 <- mean8 + se8*1.645
lower90_8 <- mean8 - se8*1.645
p8 <- summary(lm(data = min, min_ben ~ treat))$coefficients[2,4]

txdata <- read.csv("./tex.csv")
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

mean5 <- summary(lm(data = tx_pos, Tx_ben ~ treatment))$coefficients[2,1]
se5 <- summary(lm(data = tx_pos, Tx_ben ~ treatment))$coefficients[2,2]
upper95_5 <- mean5 + se5*1.96
lower95_5 <- mean5 - se5*1.96
upper90_5 <- mean5 + se5*1.645
lower90_5 <- mean5 - se5*1.645
p5 <- summary(lm(data = tx_pos, Tx_ben ~ treatment))$coefficients[2,4]

mean6 <- summary(lm(data = tx_neg, Tx_ben ~ treatment))$coefficients[2,1]
se6 <- summary(lm(data = tx_neg, Tx_ben ~ treatment))$coefficients[2,2]
upper95_6 <- mean6 + se6*1.96
lower95_6 <- mean6 - se6*1.96
upper90_6 <- mean6 + se6*1.645
lower90_6 <- mean6 - se6*1.645
p6 <- summary(lm(data = tx_neg, Tx_ben ~ treatment))$coefficients[2,4]

min <- read.csv("./min.csv")

min <- min %>% filter(treat < 3) %>% mutate(treat = treat - 1)

mean11 <- summary(lm(data = min, personal_ben ~ treat))$coefficients[2,1]
se11 <- summary(lm(data = min, personal_ben ~ treat))$coefficients[2,2]
upper95_11 <- mean11 + se11*1.96
lower95_11 <- mean11 - se11*1.96
upper90_11 <- mean11 + se11*1.645
lower90_11 <- mean11 - se11*1.645
p11 <- summary(lm(data = min, personal_ben ~ treat))$coefficients[2,4]

min <- read.csv("./min.csv")

min <- min %>% filter(treat > 2) %>% mutate(treat = treat - 3)

mean12 <- summary(lm(data = min, personal_ben ~ treat))$coefficients[2,1]
se12 <- summary(lm(data = min, personal_ben ~ treat))$coefficients[2,2]
upper95_12 <- mean12 + se12*1.96
lower95_12 <- mean12 - se12*1.96
upper90_12 <- mean12 + se12*1.645
lower90_12 <- mean12 - se12*1.645
p12 <- summary(lm(data = min, personal_ben ~ treat))$coefficients[2,4]

txdata <- read.csv("./tex.csv")
tx_pos <- txdata %>% filter(treatment < 2)
tx_neg <- txdata %>% filter(treatment > 1) %>% mutate(treatment = treatment - 2)

mean9 <- summary(lm(data = tx_pos, personal_ben ~ treatment))$coefficients[2,1]
se9 <- summary(lm(data = tx_pos, personal_ben ~ treatment))$coefficients[2,2]
upper95_9 <- mean9 + se9*1.96
lower95_9 <- mean9 - se9*1.96
upper90_9 <- mean9 + se9*1.645
lower90_9 <- mean9 - se9*1.645
p9 <- summary(lm(data = tx_pos, personal_ben ~ treatment))$coefficients[2,4]

mean10 <- summary(lm(data = tx_neg, personal_ben ~ treatment))$coefficients[2,1]
se10 <- summary(lm(data = tx_neg, personal_ben ~ treatment))$coefficients[2,2]
upper95_10 <- mean10 + se10*1.96
lower95_10 <- mean10 - se10*1.96
upper90_10 <- mean10 + se10*1.645
lower90_10 <- mean10 - se10*1.645
p10 <- summary(lm(data = tx_neg, personal_ben ~ treatment))$coefficients[2,4]

library(ggplot2)
library(ggbreak)

# insert the values into a df
plot_data <- data.frame(experiment = c("n","m","l","k","j","i","h","g","f","e","d","c","b","a"), State = c("Texas", "Texas", "Minnesota", "Minnesota",NA, "Texas", "Texas", "Minnesota", "Minnesota",NA, "Texas", "Minnesota", "Minnesota", "Minnesota"), means = c(mean1,mean3,mean2,mean4,NA,mean5,mean7,mean6,mean8,NA,mean9,mean11,mean10,mean12), pvalues= c(p1,p3,p2,p4,NA,p5,p7,p6,p8,NA,p9,p11,p10,p12), lower95 = c(lower95_1, lower95_3, lower95_2, lower95_4,NA,lower95_5, lower95_7, lower95_6, lower95_8,NA,lower95_9, lower95_11, lower95_10, lower95_12), upper95 = c(upper95_1, upper95_3, upper95_2, upper95_4,NA,upper95_5, upper95_7, upper95_6, upper95_8,NA,upper95_9, upper95_11, upper95_10, upper95_12), lower90 = c(lower90_1, lower90_3, lower90_2, lower90_4,NA,lower90_5, lower90_7, lower90_6, lower90_8,NA,lower90_9, lower90_11, lower90_10, lower90_12), upper90 = c(upper90_1, upper90_3, upper90_2, upper90_4,NA,upper90_5, upper90_7, upper90_6, upper90_8,NA,upper90_9, upper90_11, upper90_10, upper90_12))

# plot the values from the df for later label editing in pptx
ggplot(plot_data)+
  geom_pointrange(aes(x=experiment, y=means, ymin=lower95, ymax=upper95, col = State), fatten = 1)+
  geom_pointrange(aes(x=experiment, y=means, ymin=lower90, ymax=upper90, col = State), fatten = 1, size = 2)+
  ylim(c(-1.1,1.1))+
  geom_hline(yintercept = 0, linetype=2, color = 'red') +
  annotate("rect", xmin = 0.5, xmax = 2.5, ymin = -1.1, ymax = 1.1,
           alpha = .3,fill = "grey") +
  annotate("rect", xmin = 5.5, xmax = 7.5, ymin = -1.1, ymax = 1.1,
           alpha = .3,fill = "grey") +
  annotate("rect", xmin = 10.5, xmax = 12.5, ymin = -1.1, ymax = 1.1,
           alpha = .3,fill = "grey") +
  annotate("rect", xmin = 0.5, xmax = 9.5, ymin = -1.1, ymax = 1.1,
           alpha = .1,fill = "pink") +
  annotate(geom = "text", x = 14.15, y = 1.02, label = c("p = 0.0169"), size = 2.75) +
  annotate(geom = "text", x = 13.15, y = 1.02, label = c("p = 0.0672"), size = 2.75) +
  annotate(geom = "text", x = 12.15, y = 1.02, label = c("p = 0.0729"), size = 2.75) +
  annotate(geom = "text", x = 11.15, y = 1.02, label = c("p = 0.0098"), size = 2.75) +
  annotate(geom = "text", x = 9.15, y = 1.02, label = c("p = 0.0028"), size = 2.75) +
  annotate(geom = "text", x = 8.15, y = 1.02, label = c("p = 0.0001"), size = 2.75) +
  annotate(geom = "text", x = 7.15, y = 1.02, label = c("p = 0.0542"), size = 2.75) +
  annotate(geom = "text", x = 6.15, y = 1.02, label = c("p = 0.0318"), size = 2.75) +
  annotate(geom = "text", x = 4.15, y = 1.02, label = c("p = 0.3167"), size = 2.75) +
  annotate(geom = "text", x = 3.15, y = 1.02, label = c("p = 0.3599"), size = 2.75) +
  annotate(geom = "text", x = 2.15, y = 1.02, label = c("p = 0.4000"), size = 2.75)+
  annotate(geom = "text", x = 1.15, y = 1.02, label = c("p = 0.5945"), size = 2.75) +
  coord_cartesian(ylim = c(35, 65), expand = FALSE, clip = "off") +
  labs(y= "", x = "") + scale_x_discrete(labels = c("Experiment 4","Experiment 2", "Experiment 3", "Experiment 1", "","Experiment 4","Experiment 2", "Experiment 3", "Experiment 1", "","Experiment 4","Experiment 2", "Experiment 3", "Experiment 1"))+
  theme_bw() + theme(legend.position = "none") + coord_flip() + 
  scale_colour_discrete(na.translate = F)

#---

# Mediation Analysis for Figure 5

library(mediation)
library(sandwich)
library(ggplot2)

mindata <- read.csv("./min.csv")
min_neg <- mindata %>% filter(treat > 2) %>% mutate(treat = treat - 3)
min_pos <- mindata %>% filter(treat < 3) %>% mutate(treat = treat - 1)
# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 5
set.seed(5)
med.fit.postex <- lm(Tx_ben ~ treatment + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw , data = tx_pos)
out.fit.postex <- lm(support ~ treatment + Tx_ben + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw, data = tx_pos)

med.out.postex <- mediate(med.fit.postex, out.fit.postex, treat = "treatment", boot = TRUE, mediator = "Tx_ben", robustSE = FALSE)
summary(med.out.postex)

# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 2
set.seed(2)
med.fit.negtex <- lm(Tx_ben ~ treatment + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw, data = tx_neg)
out.fit.negtex <- lm(support ~ treatment + Tx_ben + Gender + Age + Education +  Income + Born.in.U.S. + Partisanship + Mexican + American.Id + Texas.Id + txknw, data = tx_neg)

med.out.negtex <- mediate(med.fit.negtex, out.fit.negtex, treat = "treatment", boot = TRUE, mediator = "Tx_ben", robustSE = FALSE)
summary(med.out.negtex)

# summary(med.out.pos)
# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 11
set.seed(11)
med.fit.posmin <- lm(min_ben ~ treat + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_pos)
out.fit.posmin <- lm(support ~ treat + min_ben + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_pos)

med.out.posmin <- mediate(med.fit.posmin, out.fit.posmin, treat = "treat", boot = TRUE, mediator = "min_ben", robustSE = FALSE)
summary(med.out.posmin)

# Random number func: runif(n=1, min=1, max=20)
# Source: https://www.statology.org/r-random-number/
# Random number: 17
set.seed(17)
med.fit.negmin <- lm(min_ben ~ treat + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_neg)
out.fit.negmin <- lm(support ~ treat + min_ben + Gender + Age + education +  income + us_born + partisanship + white + American.Id + Minnesota.Id + minkn, data = min_neg)

med.out.negmin <- mediate(med.fit.negmin, out.fit.negmin, treat = "treat", boot = TRUE, mediator = "min_ben", robustSE = FALSE)
summary(med.out.negmin)

